﻿using System;
using LatoolNet;

namespace NumericalSolution {
  class Poisson2DFVDirect : Poisson2DFV {

    private Matrix m_u;
    private Matrix m_b;

    public override void Prepare() {

      int rownum = p_numberOfControlVolumes * p_numberOfControlVolumes;
      int colnum = p_numberOfControlVolumes * p_numberOfControlVolumes;
      int bandWidth = 2 * p_numberOfControlVolumes + 1;
      m_u = new Matrix(rownum, colnum, bandWidth);
      m_b = new Matrix(rownum, 1);
    
    }

    public override ROResult2D Solve() {

      for (int i = 0; i < p_numberOfControlVolumes; i++) {
        for (int j = 0; j < p_numberOfControlVolumes; j++) {

          double coef_w = Double.NaN;
          double coef_e = Double.NaN;
          double coef_n = Double.NaN;
          double coef_s = Double.NaN;
          double coef_p = 0.0;

          if (i == 0) {
            coef_p += 2 * p_deltay / p_deltax;
            coef_w = 0.0;
          } else {
            coef_p += p_deltay / p_deltax;
            coef_w = p_deltay / p_deltax;
          }

          if (i == p_numberOfControlVolumes - 1) {
            coef_p += 2 * p_deltay / p_deltax;
            coef_e = 0.0;
          } else {
            coef_p += p_deltay / p_deltax;
            coef_e = p_deltay / p_deltax;
          }

          if (j == 0) {
            coef_p += 2 * p_deltax / p_deltay;
            coef_s = 0.0;
          } else {
            coef_p += p_deltax / p_deltay;
            coef_s = p_deltax / p_deltay;
          }

          if (j == p_numberOfControlVolumes - 1) {
            coef_p += 2 * p_deltax / p_deltay;
            coef_n = 0.0;
          } else {
            coef_p += p_deltax / p_deltay;
            coef_n = p_deltax / p_deltay;
          }

          m_u[index(i, j), index(i, j)] = coef_p;

          if (i != 0) {
            m_u[index(i, j), index(i - 1, j)] = -coef_w;
          }
          if (i != p_numberOfControlVolumes - 1) {
            m_u[index(i, j), index(i + 1, j)] = -coef_e;
          }
          if (j != p_numberOfControlVolumes - 1) {
            m_u[index(i, j), index(i, j + 1)] = -coef_n;
          }
          if (j != 0) {
            m_u[index(i, j), index(i, j - 1)] = -coef_s;
          }
          m_b[index(i, j), 0] = f(posx(i), posy(j)) * p_deltax * p_deltay;

        }
      }
      //LUFactorization.Solve(m_u, m_b);
      LinearEquation.Solve(m_u, m_b);

      return new ROResult2D(m_b, p_numberOfControlVolumes);
    }

  }
}
